library(dplyr)
library(readr)
library(lubridate)
library(zoo)
library(sandwich)

## interest rate -- fed funds rate
data <- read_csv("data/FEDFUNDS.csv") %>%
    rename(date = DATE, i = FEDFUNDS) %>%
    mutate(date = as.Date(date, format = "%m/%d/%Y"),
           ym = year(date)*100+month(date))
## lagged interest rate
data <- mutate(data, l.i = dplyr::lag(i, 3))

## output gap
gap <- read_csv("data/output_gap.csv") %>%
    rename(date = DATE, gap = GDPC1_GDPPOT)
data <- left_join(data, gap) %>%
    mutate(gap = na.approx(gap, na.rm=FALSE),
           gap = na.locf(gap))

## inflation
inf <- read_csv("data/CPIAUCSL.csv") %>% rename(date = DATE, P = CPIAUCSL) %>%
    mutate(p = log(P),
           inf =  100*(p - dplyr::lag(p, 12))) %>%
    select(date, inf)
data <- left_join(data, inf) %>%
    na.omit

## rolling regression
t0 <- which(data$ym == 198501)
## window length
m <- 7*12   # seven years
data <- data %>%
    mutate(gamma = NA, beta = NA, rho = NA, se_gamma = NA, se_beta = NA, se_rho = NA)
for (t in t0:nrow(data)) {
    ## rolling-window least squares
    ind <- max(1, t-m+1):t
    mod <- lm(i ~ l.i + gap + inf, data[ind,])
    ## save coefficients
    data$gamma[t] <- mod$coef["gap"]
    data$beta[t] <- mod$coef["inf"]
    data$rho[t] <- mod$coef["l.i"]
    SEs <- sqrt(diag(NeweyWest(mod, prewhite=FALSE, lag=12)))
    data$se_gamma[t] <- SEs["gap"]
    data$se_beta[t] <- SEs["inf"]
    data$se_rho[t] <- SEs["l.i"]
}
data$gamma_lb <- data$gamma - 2*data$se_gamma
data$gamma_ub <- data$gamma + 2*data$se_gamma
data$beta_lb <- data$beta - 2*data$se_beta
data$beta_ub <- data$beta + 2*data$se_beta
data$rho_lb <- data$rho - 2*data$se_rho
data$rho_ub <- data$rho + 2*data$se_rho

cat("First window:", format(range(data$date[(t0-m+1):t0])), "\n")
cat("Last window:", format(range(data$date[ind])), "\n")

data <- data %>% filter(year(date) >= 1985)

data <- data %>%
    select(date, beta, gamma, rho,
           beta_lb, beta_ub, gamma_lb, gamma_ub, rho_lb, rho_ub) %>%
    rename(beta_actual = beta,
           beta_actual_lb = beta_lb,
           beta_actual_ub = beta_ub,
           gamma_actual = gamma,
           gamma_actual_lb = gamma_lb,
           gamma_actual_ub = gamma_ub,
           rho_actual = rho,
           rho_actual_lb = rho_lb,
           rho_actual_ub = rho_ub)
write.csv(data, file="output/actual_rule_inertial.csv", row.names=FALSE, quote=FALSE)
save(data, file="output/actual_rule_inertial.RData")
